home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / x / gui / xfract.lha / xfract / gen_mntndata.c < prev    next >
Encoding:
C/C++ Source or Header  |  1990-06-27  |  4.3 KB  |  165 lines
  1. /****************************************************************************
  2. /* FILE:    gen_mntndata.c
  3. /* AUTHOR:    Paul Sharpe @ DEC, OSCR_Europe, reading, England.
  4. /* DATE:    October 11, 1989.
  5. /* INSPIRATION:    'The Science Of factal Images'
  6. /*
  7. /*   Copyright (c) Digital Equipment Corporation 1990  All rights reserved.
  8. /*   Copyright is claimed in the computer program and user interface thereof.
  9. /*
  10. /*   Digital Equipment Corporation cannot accept any responsibility for
  11. /*   use, misuse, or abuse of this software.
  12. /*
  13. /****************************************************************************/
  14.  
  15. #include <stdio.h>
  16. #include <math.h>
  17.  
  18. #include "xpt.h"
  19.  
  20. #define F3(delta,x0,x1,x2)      ((x0+x1+x2)/3 + delta*Gauss())
  21. #define F4(delta,x0,x1,x2,x3)   ((x0+x1+x2+x3)/4 + delta*Gauss())
  22.  
  23. int        level;
  24. char        addition = 1;
  25. double        sigma, H;
  26. static int    Nrand, Arand;
  27. static double    GaussFac, GaussAdd;
  28. double        Gauss();
  29. static double    numer, denom; /* Replaces use of GaussFac for precision... */
  30.  
  31. gen_mntndata(elevs, square, top)
  32. double    **elevs, top;
  33. int    square;
  34. {
  35.     MidPointFM2D(elevs,level,sigma,H,addition,(int)time((long *)0));
  36.     scale(elevs, square, top);
  37. }
  38.  
  39. scale(elevs, square, top)
  40. double    **elevs;
  41. int    square;
  42. double    top;
  43. {
  44. register int    x,y;
  45. double        max = 0.0;
  46.  
  47.     for (x=0; x<square; x++)
  48.     for (y=0; y<square; y++)
  49.         if (elevs[x][y] > max)
  50.         max = elevs[x][y];
  51.     DEBUG(("Scaling: %f -> %f.\n",max,top));
  52.     if (max > 0.0)            /* Just in case summat goes wrong... */
  53.     for (x=0; x<square; x++)
  54.         for (y=0; y<square; y++)
  55.         elevs[x][y] = elevs[x][y]*top/max;
  56. }
  57.  
  58. /*==========================================================================
  59.  * Routines 'interpreted' from the code in 'The Science Of Fractal Images.'
  60.  *==========================================================================*/
  61. InitGauss(seed)
  62. int    seed;
  63. {
  64.     Nrand = 4;
  65.     Arand = (1<<31)-1;        /* random() gives results in [0,Arand] */
  66.  
  67.     GaussAdd = (double)sqrt((double)(3.0*Nrand));
  68. /* NB: GaussFac is no longer used: it loses too much precision!
  69.  *   GaussFac = (double)(2.0 * GaussAdd / ((double)Nrand * (double)Arand));
  70.  */
  71. /* NB: We have to do this, as we seem to loose too much precision in the method
  72.  *     laid out in 'The Scienc Of Fractal Images.'
  73.  */
  74.     numer = GaussAdd + GaussAdd;
  75.     denom = (double)((double)Nrand * (double)Arand);
  76.  
  77.     srandom(seed);
  78. }
  79.  
  80. double
  81. Gauss()
  82. {
  83. register int    i;
  84. register double    sum = 0;
  85. extern long    random();
  86.  
  87. /* See the above comment in InitGauss()... */
  88.     for (i=0; i<Nrand; i++)
  89.     sum += (double)random();
  90.     sum = numer*(sum/denom);
  91.     return(sum - GaussAdd);
  92. }
  93.  
  94. MidPointFM2D(X,maxlevel,sigma,H,addition,seed)
  95. double    **X;
  96. int    maxlevel;
  97. double    sigma;
  98. double    H;
  99. char    addition;
  100. int    seed;
  101. {
  102. register int    i, N, x,y, D,d;
  103. double        delta;
  104.  
  105.     InitGauss(seed);
  106.     N = 1<<maxlevel;
  107.     DEBUG(("MidPointFM2D: (Dim:%dx%d) (Sigma:%f) (H:%f).\n",N,N,sigma,H));
  108.  
  109.     delta = sigma;
  110.     X[0][0] = delta * Gauss();
  111.     X[0][N] = delta * Gauss();
  112.     X[N][0] = delta * Gauss();
  113.     X[N][N] = delta * Gauss();
  114.  
  115.     D = N;
  116.     d = N/2;
  117.  
  118.     for (i=0; i<maxlevel; i++) {
  119.  
  120. /* Going from Grid type I to Grid type II */
  121.     delta = delta * pow((double)0.5,(double)(0.5*H));
  122.     for (x=d; x<=N-d; x+=D)        /* Interpolate and offset points. */
  123.         for (y=d; y<=N-d; y+=D)
  124.         X[x][y] = F4(delta,X[x+d][y+d],X[x+d][y-d],X[x-d][y+d],X[x-d][y-d]);
  125.  
  126. /* Displace other points also if needed. */
  127.     if (addition == 1)
  128.         for (x=0; x<=N; x+=D)
  129.         for (y=0; y<=N; y+=D)
  130.             X[x][y] = X[x][y] + delta*Gauss();
  131.  
  132. /* Going from Grid type II to Grid type I */
  133.     delta = delta * pow(0.5,0.5*H);
  134. /* Interpolate and offset boundary grid points. */
  135.     for (x=d; x<=N-d; x+=D) {
  136.         X[x][0] = F3(delta,X[x+d][0],X[x-d][0],X[x][d]);
  137.         X[x][N] = F3(delta,X[x+d][N],X[x-d][N],X[x][N-d]);
  138.         X[0][x] = F3(delta,X[0][x+d],X[0][x-d],X[d][x]);
  139.         X[N][x] = F3(delta,X[N][x+d],X[N][x-d],X[N-d][x]);
  140.     }
  141.  
  142. /* Interpolate and offset interior grid points. */
  143.     for (x=d; x<=N-d; x+=D)
  144.         for (y=D; y<=N-d; y+=D)
  145.         X[x][y] = F4(delta,X[x][y+d],X[x][y-d],X[x+d][y],X[x-d][y]);
  146.  
  147.     for (x=D; x<=N-d; x+=D)
  148.         for (y=d; y<=N-d; y+=D)
  149.         X[x][y] = F4(delta,X[x][y+d],X[x][y-d],X[x+d][y],X[x-d][y]);
  150.  
  151. /* Displace other points also if needed. */
  152.     if (addition == 1) {
  153.         for (x=0; x<=N; x+=D)
  154.         for (y=0; y<=N; y+=D)
  155.             X[x][y] = X[x][y] + delta*Gauss();
  156.         for (x=d; x<=N-d; x+=D)
  157.         for (y=d; y<=N-d; y+=D)
  158.             X[x][y] = X[x][y] + delta*Gauss();
  159.     }
  160.  
  161.     D = D/2;
  162.     d = d/2;
  163.     }
  164. }
  165.